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Abstract 

A simple and explicit technique for the numerical solution of the two-particle, time- 
dependent Schrodinger equation is assembled and tested. The technique can handle 
interparticle potentials that are arbitrary functions of the coordinates of each parti- 
cle, arbitrary initial and boundary conditions, and multi-dimensional equations. Plots 
and animations are given here and on the World Wide Web of the scattering of two 
wavepackets in one dimension 

1 Introduction 

Rather than showing the time dependence of two particles interacting with each other, 
quantum mechanics textbooks often present a time-independent view of a single particle 
interacting with an external potential. In part, this makes the physics clearer, and in part, 
this reflects the difficulty of solving the time-independent two-particle Schrodinger equation 
for the motion of wavepackets. In the classic quantum mechanics text by Schiff Q, examples 
of realistic quantum scattering, such as that in Fig. |^, are produced by computer simulations 
of wave packets colliding with square potential barriers and wells. Generations of students 
have carried memories of these images (or of the film loops containing these frames [^) as 
to what realistic quantum scattering looks like. 

While Fig. | is a good visualization of a quantum scattering processes, we wish to 
extend simulations of realistic quantum interactions to include particle-particle scattering 
when both particles are represented by wavepackets. Although more complicated, this, 
presumably, is closer to nature and may illustrate some physics not usually found in quantum 
mechanics textbooks. In addition, our extension goes beyond the treatment found in most 
computational physics texts which concentrate on one-particle wavepackets [^, ^, ^, or 
highly restricted forms of two-particle wavepackets [^. 

The simulations of the time-dependent Schrodinger equation shown by Schiff were based 
on the 1967 finite-difference algorithms developed by Goldberg et al. Those simulations, 
while revealing, had problems with stability and probability conservation. A decade later, 
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Figure 1: A time sequence of a Gaussian wavepacket scattering from a square barrier as 
taken from the textbook by Schiff. The mean energy equals the barrier height. 



Cakmak and Askar solved the stability problem by using a better approximation for the 
time derivative. After yet another decade, Visscher solved the probability conservation 
problem by solving for the real and imaginary parts of the wave function at slightly different 
("staggered") times. 

In this paper we combine the advances of the last 20 years and extend them to the 
numerical solution of the two •particle — in contrast to the one "particle — time-dependent 
Schrodinger equation. Other than being independent of spin, no assumptions are made re- 
garding the functional form of the interaction or initial conditions, and, in particular, there 
is no requirement of separation into relative and center-of-mass variables]^]. The method is 
simple, explicit, robust, easy to modify, memory preserving, and may have research applica- 
tions. However, high precision does require small time and space steps, and, consequently, 
long running times. A similar approach for the time-dependent one-particle Schrodinger 
equation in a two-dimensional space has also been studied [^. 



2 Two-Particle Schrodinger Equation 

We solve the two-particle time-dependent Schrodinger equation 
d 

i-Q^'4'{xi,X2,t) = H'llj{xi,X2,t), (1) 

1 52 1 92 

2midx'\ 2m2 ~'~ ^ ^ 

where, for simplicity, we assume a one-dimensional space and set fi = \. Here H is the 
Hamiltonian operator and rrii and Xi are the mass and position of particle i = 1, 2. Knowl- 
edge of the two-particle wave function ■0(2^1, 2^25 permits the calculation of the probability 
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Figure 2: Six frames from an animation of the two-particle density p{xi,X2,t) as a function 
of the particle positions xi and X2, for a repulsive m-lOm collision in which the mean kinetic 
energy equals twice the barrier height. The numbers in the left hand corners give the times 
in units of lOOAt. 



density for particle 1 being at xi and particle 2 being at X2 at time t: 

p{xi,X2,t) = \ip{xi,X2,t)\'^ . (3) 

The fact that particles 1 and 2 must be located someplace in space leads to the normalization 
constraint on the wave function: 

/ dxidx2\i}ixi,X2,t)\ =1. (4) 

-oo J —oo 

The description of a single particle within a multi-particle system by a single-particle 
wave function is an approximation unless the system is uncorrelated (in which case the total 
wave function can be written in product form). However, it is possible to deduce meaningful 
one-particle densities from the two-particle density by integrating over the other particle: 

/ + 00 
dxj p{xi,X2,t), {i^j = 1,2). (5) 

Here we use a subscript on the single-particle density pi to distinguish it from the two- 
particle density p. Of course, the true solution is ip{xi,X2,t), but we find it hard to see the 
physics in a three- variable complex function, and so, often, view pi{xi,t) and P2{x2,t) as 
two separate wavepackets colliding. 
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If particles 1 and 2 are identical, then their total wave function should be symmetric 
or antisymmetric under interchange of the particles. We impose this condition on our 
numerical solution ip(xi,X2), by forming the combinations 

'4^'{xi,X2) = [^/^(xi,a;2) ± V'(x2,xi)] (6) 

2p{xi,X2) = \tp{xi,X2)\'^ + \ip{x2,xi)\'^ ±2Re [■ip*{xi,X2)ip{x2,xi)] . (7) 
The cross term in places an additional correlation into the wavepackets. 



3 Numerical Method 

We solve the two-particle Schrodinger equation (jl]) via a finite difference method that con- 
verts the partial differential equation into a set of simultaneous, algebraic equations. First, 
we evaluate the dependent variable -0 on a grid of discrete values for the independent vari- 
ables §: 

xp{xi,X2,t) = ip{xi = lAxi,X2 = mAx2,t = nAt) = (8) 

where I, m, and n are integers. The space part of the algorithm is based on Taylor expansions 
of 'ip{xi, X2,t) in both the xi and X2 variables up to 0{Ax^); for example, 

d'^lp Ipixi + Axi,X2) -2ll}{xi,X2) + Tp{xi - Axi,X2) 2n /n\ 

TTT TT^ + 0{Ax^). (9 

oxi Axf 

In discrete notation, the RHS of the Schrodinger equation (|^) now becomes: 

= 2^^, 2^^, + 

Next, we express the time derivative in (||) in terms of finite time differences by taking the 
formal solution to the time-dependent Schrodinger equation and making a forward-difference 
approximation for time evolution operator: 

= e-^^^^V'zV =^ il-^AtH)i;r,m■ (H) 

Although simple, this approximation scheme is unstable since the term multiplying tp has 
eigenvalue (1 — iEAt) and modulus Vl + E^At"^, and this means the modulus of the wave 
function increases with each time step The improvement introduced by Askar and 
Cakmak [Q] is a central difference algorithm also based on the formal solution ([ll|): 



= 




l,m 




n+l ^ 




l,in 














[mi 



\{± + ±)AX + AxViA4:r,m (13) 



' ' 772-2 

where we have assumed Axi = Ax2 and formed the ratio A = At/Ax'^. 

Equation (0) is an explicit solution in which the wave function at only two past time 
values must be stored simultaneously in memory to determine all future times by continued 
iteration. In contrast, an implicit solution determines the wave function for all future times 
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in just one step, yet this one step requires the solution of simultaneous algebraic equations 
involving all space and time values. Accordingly, an implicit solution requires the inversion 
of exceedingly large (~ 10^'^ x 10^'') matrices. 

While the explicit method ( [l3|) produces a solution which is stable and second-order 
accurate in time, in practice, it does not conserve probability well. Visscher|^] has deduced 
an improvement which takes advantage of the extra degree of freedom provided by the 
complexity of the wave function to preserve probability better. If we separate the wave 
function into real and imaginary parts. 



n+l _ ,,n+l 



Lm 
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the algorithm (13) separates into the pair of coupled equations: 
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(14) 



(15) 



(16) 



Visscher's advance evaluates the real and imaginary parts of the wave function at slightly 
different (staggered) times. 



Km, <m] = [Re'^{x,t), ImV'(x,t + ^At)], 



(17) 



and uses a definition for probability density that differs for integer and half-integer time 
steps. 



pix,t) 
At^ 



\Reilj[x,t)\ +lmil;{x,t + —)lmip{x,t —), 



p{x,t + —) = Re'ilj{x,t + At)Reilj{x,t) + 



At 2 

Imi/;(x, t H — —) 



(18) 
(19) 



These definitions reduce to the standard one for infinitesimal At, and provide an algebraic 
cancellation of errors so that probability is conserved. 



4 Simulations 

We assume that the particle-particle potential is central and depends only on the relative 
distance between particles 1 and 2 (the method can handle any xi and X2 functional de- 
pendences). We have investigated a "soft" potential with a Gaussian dependence, and a 
"hard" one with a square- well dependence, both with range a and depth Vq: 

V{xuX2) = /^oexp[-^^;^] (Gaussian) _ ^20) 
[Vo 9{a — \xi — X2\) (Square) 
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Table 1: Table 1, Parameters for the antisymmetrized, m-m collision with an attractive 
square well potential. 
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Figure 3: A time sequence of two Gaussian single-particle wavepackets scattering from each 
other under the influence of a square barrier. The mean energy equals twice the barrier 
height. The dashed curve describes a particle of mass m and the solid curve one of mass 
10m. The number in the upper left-hand corner of each frame is the time in units of lOOAi, 
and the edges of the frames correspond to the walls of the box. 
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4.1 Initial and Boundary Conditions 

We model a scattering experiment in which particle 1, initially at with momentum ki, 
collides with particle 2, initially far away at X2 with momentum k2, by assuming a product 
of independent wavepackets for particles 1 and 2: 



V'(xi,X2,i = 0) 



{xi 



X 



0^2 



4a2 



{X2 



4C72 



(21) 



Because of these Gaussian factors, is not an eigenstate of the particle i momentum 
operators —id/dxi, but instead contains a spread of momenta about the mean, initial 
momenta ki and ^2- If the wavepacket is made very broad (a 00), we would obtain 
momentum eigenstates. Note, that while the Schrodinger equation may separate into one 
equation in the relative coordinate x and another in the center-of-mass coordinate X, the 



initial condition (21), or more general ones, cannot be written as separate conditions on x 
and X. Accordingly, a solution of the equation in each particle coordinate is required |^]. 

We start the staggered-time algorithm with the real part the wave function (21) at t = 
and the imaginary part at t = At/ 2. The initial imaginary part follows by assuming that 
At/2 is small enough, and a large enough, for the initial time dependence of the wavepacket 
to be that of the plane wave parts: 



lm'ip{xi,X2,t = — ) 



sm 



X exp 




(22) 



In a scattering experiment, the projectile enters from infinity and the scattered particles 
are observed at infinity. We model that by solving our partial differential equation within 
a box of side L (ideally) much larger than both the range of the potential and the width of 
the initial wavepacket. This leads to the boundary conditions 



^'(O, X2,t) = 'ip{xi,0, t) = ip{L, X2,t) = ip{xi,L, t) = 0. 



(23) 



The largeness of the box minimizes the effects of the boundary conditions during the collision 
of the wavepackets, although at large times there will be interesting, yet artificial, collisions 
with the box. 

Some typical parameters used in our tests are given in Table |l| (the code with sample 
files are available on the on Web Q). Our space step size Ax = 0.001 is 1/1, 400th of 
the size of the box L, and l/70th of the size {ypla ~ 0.07) of the wavepacket. Our time 
step At = 2.5 X 10"^ is 1/20, 000th of the total time T, and 1/2, 000th of a typical time 
for the wavepacket [27r/(/cf/2mi) ~ 5 x 10^^]. In all cases, the potential and wavepacket 
parameters are chosen to be similar to those used in the one-particle studies by Goldberg 
et al. The time and space step sizes were determined by trial and error until values were 
found which provided stability and precision (too large a Ax leads to spurious ripples 
during interactions). In general, stability is obtained by making At small enough with 
simultaneous changes in At and Ax made to keep A = At/Ax^ constant. Total probability, 
as determined by a double Simpson's-rule integration of (Q), is typically conserved to 13 
decimal places, impressively close to machine precision. In contrast, the mean energy, for 
which we do not use a definition optimized to staggered times, is conserved only to 3 places. 
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Figure 4: Same as Fig. ^ except now the potential is attractive with the mean energy equal 
to half the depth. 

4.2 Barrier-Like Collisions 

We solve our problem in the center-of-momentum system by taking k2 = —ki (particle 
1 moving to larger x values and particle 2 to smaller x). Our first simulations and Web 
animations emulate the one-particle collisions with barriers and wells studied by Goldberg 
et al. and presented by Schiff. We make particle 2 ten times heavier than particle 1, so that 
particle 2's initial wavepacket moves at 1/lOth the speed of particle I's, and so looks like 
a barrier. Although we shall describe several scattering events, the animations available on 
the Web speak for themselves, and we recommend their viewing. 

In Fig. I we show six frames from an animation of the two-particle density p{xi,X2,t) 
as a simultaneous function of the particle positions xi and X2. In Fig. ^ we show, for 
this same collision, the single-particle densities pi{x = xi,t) and p2{x = X2,t) extracted 
from p{xi,X2,t) by integrating out the dependence on the other particle via (^). Since the 
mean energy equals twice the maximum height of the potential barrier, we expect complete 
penetration of the packets, and indeed, at time 18 we see that the wavepackets have large 
overlap, with the repulsive interaction "squeezing" particle 2 (it gets narrower and taller). 
During times 22-40 we see part of wavepacket 1 reflecting off wavepacket 2 and then moving 
back to smaller x (the left). From times 26-55 we also see that a major part of wavepacket 
1 gets "trapped" inside of wavepacket 2 and then leaks out rather slowly. 

We see that for times 1-26, the X2 position of the peak of p{xi,X2,t) in Fig. || changes 
very little with time, which is to be expected since particle 2 is heavy. In contrast, the 
xi dependence in p{xi,X2,t) gets broader with time, develops into two peaks at time 26, 
separates into two distinct parts by time 36, and then, at time 86 after reflecting off the 
walls, returns to particle 2's position. We also notice in both these figures, that at time 
40 and thereafter, particle 2 (our "barrier" ) fissions and there are two peaks in the X2 
dependence. 

As this comparison of Figures ^ and |3| demonstrates, it seems easier to understand the 
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m-m, Repulsive Vsquare 
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Figure 5: Same as Fig. ^, except now for a repulsive m-m collision in which the mean 
energy equals one quarter of the barrier's height. 

physics by superimposing two single-particle densities (thereby discarding information on 
correlations) than by examining the two-particle density. Accordingly, the figures we show 
hence, and the majority of the animations on the Web, are of single-particle densities. 

Fig. ^ is similar to the behavior present in Schiff 's one-particle simulation. Fig. ||, but 
without ripples during the collision. Those ripples are caused by interference between 
scattered and incident wave, and even though we have a square barrier potential acting 
between the particles, neither particle "feels" the discontinuity of the sharp potential edge 
at any one time. However, there are ripples when our wavepackets hit the walls, as seen at 
times 55 and 150. (There is also a ripple at time 36 arising from interference with the small, 
reflected part of the left edge of I's initial wavepacket. This artifact of having particle 1 so 
near the left edge of the bounding box can be seen reflecting off the left wall at time 18.) 

Something new in Fig. ^, that is not in Schiff, is the delayed "fission" of the heavier 
particle's wavepacket after time 40 due to repulsion from the reflected and transmitted parts 
of wavepacket 1. In addition, at time 86 we see that the reflected and transmitted parts of 
1 have reconstituted themselves into a single but broadened wavepacket, that at time 150 
is again being reflected from the left wall. 

In Fig. I we see another m-lOm collision. This time there is an attractive interaction 
between the particles and again the mean energy equals half the well depth. Even though 
the kinetic energy is low, the interaction is attractive and so particle 1 passes through 
particle 2. However, some of wavepacket 1 is reflected back to the left after the collision, 
and, as we see at time 55, the wavepacket for the heavy particle 2 fissions as a consequence 
of its attraction to the two parts of wavepacket 1. 

Although we do not show them here, on the Web we also display movies of collisions 
corresponding to a Gaussian potential acting between the particles. These are much softer 
collisions and have behaviors similar to classical particles bouncing off each other, with 
squeezing and broadening of the wavepackets, but little breakup or capture. 
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m-IOm, Attractive Vsquare 
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Figure 6: Same as Fig. ^ except now for an attractive m-m collision in which the mean 
energy equals one quarter of the well's depth. 
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Figure 7: Same as Fig. ^ except now for an attractive m-m collision in which the mean 
energy equals one quarter of the well's depth, and for which the wavefunction has been 
antisymmetrized . 
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Symmetrized m-m, Attractive Vsquare 

KE = -V/4 




Figure 8: Same as Fig. except now for an attractive m-m collision in which the mean en- 
ergy equals one quarter of the depth, and for which the wavefunction has been symmetrized. 



4.3 m— m Collisions 

In Fig. |5| we show nine frames from the movie of a repulsive m-m collision in which the 
mean kinetic energy equals one quarter of the barrier height. The initial packets are seen to 
slow down as they approach each other, with their mutual repulsion narrowing and raising 
the packets up until the time (50) when they begin to bounce back. The wavepackets at still 
later times are seen to retain their shape, with a progressive broadening until they collide 
with the walls and break up. As shown on the Web, when the mean energy is raised there 
will be both transmitted and reflected wave, already seen in Fig. ^ for an m-lOm collision. 

In Fig. I we show nine frames from the movie of an attractive m-m collision in which 
the mean energy equals one quarter of the well depth. The initial packets now speed up 
as they approach each other, and at time 60 the centers have already passed through each 
other. After that, a transmitted and reflected wave for each packet is seen to develop (times 
66-78). Finally, each packet appears to capture or "pickup" a part of the other packet and 
move off with it (times 110-180). 

In Fig. ^ we repeat the collision of Fig. ^, only now for a wave function that has been 
antisymmetrized according to (0). The anitsymmetrization is seen to introduce an effective 
repulsion into what is otherwise an attraction (compare the two figures for times 60-66). 
Again, some capture of the other wavepacket is noted from times 94 on, only now the internal 
captured wavepacket retains its Gaussian-like shape, apparently the result of decreased 
interference. 

Finally, in Fig. ^ we repeat the collisions of Figures ^ and 0, only now for a wave function 
that has been symmetrized according to (|^. The symmetrization is seen to introduce an 
effective added attraction (compare the three figures for time 60 which shows the greatest 
penetration for the symmetrized case). While there is still capture of the other wavepacket, 
the movie gives the clear impression that the wavepackets interchange with each other as a 
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consequence of the symmetrization. 



5 Summary and Conclusions 

We have assembled and tested a general technique for the numerical solution of the two- 
particle, time-dependent Schrodinger equation. Because the technique is general, applica- 
tion to two or three dimensions and for other potentials and initial conditions should be 
straightforward. For example, further studies may want to investigate initial conditions 
corresponding to bound particles interacting with a surface, or the formation of a molecule 
near a surface. 

The Goldberg-Schiff 's image (Fig. ^ of a wavepacket interacting with a potential barrier 
is still a valuable model for understanding the physics occuring during a particle's collision. 
Here we have extended the level of realism to what a collision between two particles looks 
like. In doing so with a simple square-well potential between the two particles, we have 
discovered that fission and particle pickup occur quite often, although the physics may be 
quite different from that in nuclei. While somewhat of a challenge to understand fully, we 
have also provided animations of the behavior of the two-particle density during collisions. 
We have placed the animations, source codes, and movie-making instructions on the Web 
with the hope that future students will also carry some of these images of the quantum 
world with them. 
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